function [idat sandbox scaninfo] = fsb_load_files(file_names,diary,datamode)

% FSB : load files into sandbox
%
% EXAMPLE:
% [idat sandbox] = fsb_load_files(file_names,0,1)
%
% INPUT:
% sandbox.intrial:  first column refers to the trial being on or off
%                   second column refers to the number of the volume within a trial
%                   third column refers to the number of the trial
%                   fourth column refers to the number of the scan
%                   fifth column refers to the number of the session.
%                   This is determined
%                   by the first five letters of the filename.
%                   If these are identical, both are assumed to belong to
%                   the same set of files
% diary:            set this to 1 to document how you proceed processing the
%                   data
% datamode:         set this to 0 to force int16 conversion
%
% OUTPUT:
% idat:             4D image data
% sandbox:          Sandbox experiment struct
% scaninfo:         Scan data
%
% CALLED BY:
% FSB.m
% fsb_realign_multiple.m
% fsb_prepare_SPM.m
% fsb_realign_SPM.m
% fsb_calc_cluster_tc.m
%
% NOTES:
% there remains still a lot of work to be done. For backwards compatibility,
% it would be nice to be also able to load tempdat.onset and so on, but this
% has not yet been implemented
% The file names should follow the convention Subjectname: 3
% Letters/numbers, e.g. B04, scanname 2 letters (e.g. FH), and then any
% combination of numbers and letters. The full scanname could thus be
% B04FH_whatever
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
% $ Revision 1.1
%~~~~~~~~~~~~~~~~~~~~~~

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if file_names are there and call uipickfiles if needed
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<3
    datamode = 1;
end

if ~exist('file_names','var')
    file_names = uipickfiles('Output','cell','REF','.mat|.hdr|.nii');
end

if isempty(file_names)
    file_names = uipickfiles('Output','cell','REF','.mat|.hdr|.nii');
end

if ~exist('file_names','var')
    disp('No file selected, aborting...')
    return
end

keepname = 1;
subjectname = 1;
filnum = size(file_names,2);
waitstr = ['Loading ' num2str(filnum) ' files'];

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check which file_type is being loaded
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if ~isnumeric(file_names);
    
    [pathstr{1}, name{1}, ext, versn] = fileparts(file_names{1});

    filnum = size(file_names,2);

    if filnum>0

        %~~~~~~~~~~~~~~~~~~~~
        % Bruker files in 2dseq format
        %~~~~~~~~~~~~~~~~~~~~
        if strcmp(strcat(name{1}),'2dseq') == 1 ;
            file_type = '2dseq';
        end

        %~~~~~~~~~~~~~~~~~~~~
        % Analyze files in .hdr/.img format
        %~~~~~~~~~~~~~~~~~~~~
        if strcmp(strcat(ext),'.hdr') == 1 || strcmp(strcat(ext),'.img')==1;
            file_type = 'analyze';
        end
        
         %~~~~~~~~~~~~~~~~~~~~
        % Analyze files in .hdr/.img format
        %~~~~~~~~~~~~~~~~~~~~
        if strcmp(strcat(ext),'.nii') == 1;
            file_type = 'Nifti';
        end

        %~~~~~~~~~~~~~~~~~~~~
        % Matlab Arrays in .mat format
        %~~~~~~~~~~~~~~~~~~~~
        if strcmp(strcat(ext),'.mat')==1;
            file_type = 'mat';
        end

        %~~~~~~~~~~~~~~~~~~~~
        % SPM.mat file
        %~~~~~~~~~~~~~~~~~~~~
        if  strcmp(strcat(name),'SPM')==1;
            file_type = 'SPM';
        end

        if exist('diary','var') && diary == 1
            %~~~~~~~~~~~~~~~~~~~~
            % Start diary
            %~~~~~~~~~~~~~~~~~~~~
            diary(name{1});
            t = now;
            disp(['Log from ' datestr(t)]);

        end

        % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Load files
        % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        switch file_type;

            case 'analyze'
                % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Load analyze images
                % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                waitstr = ['Loading ' num2str(filnum) ' files'];
                h = waitbar(0,waitstr);

                nidat = avw_img_read_mt(strcat(file_names{1}));
                [a b c] = size(nidat.img);
                
                datatype = sub_get_data_range(nidat.img,datamode);
                idat = eval([datatype '(zeros(a,b,c,filnum))']);
                idat(:,:,:,1) = eval([datatype '(nidat.img)']);
                
                clear nidat;

                avw = avw_hdr_read(strcat(file_names{1}));
                hdr = avw.hdr;
                sandbox.hdr = hdr;
                scaninfo.hdr = hdr;
                sandbox.hemodynamics = [];
                sandbox.onset = [];
                sandbox.intrial = [];
                waitstr = (['reading file 1 of ' num2str(filnum)]);

                for x = 2:filnum
                    waitbar(x/filnum,h,waitstr)
                    nidat = avw_img_read_mt(strcat(file_names{x}));
                    waitstr = (['reading file ' num2str(x) ' of ' num2str(filnum)]);
                    disp (['Loading file ' file_names{x}]);

                    
                    nidat = eval([datatype '(nidat.img)']);
                    idat(:,:,:,x) = nidat;
                    clear nidat;
                end
                
            case 'Nifti'
                % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Load Nifti images
                % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                waitstr = ['Loading ' num2str(filnum) ' files'];
                h = waitbar(0,waitstr);

                try
                    nidat = load_nii(strcat(file_names{1}));
                catch
                    disp('Spatial transform detected in .nii image, please wait for reslicing');
                    old_fn = file_names{1};
                    new_fn = [pathstr{1} '\transform_' name{1} ext];
                    reslice_nii(old_fn, new_fn);
                    nidat = load_nii(new_fn);
                    disp('...done');
                end
                            
                [a b c] = size(nidat.img);

                datatype = sub_get_data_range(nidat.img,datamode);
                try
                    idat = nidat.img;
                catch

                    idat = eval([datatype '(zeros(a,b,c,filnum))']);
                    idat(:,:,:,1) = eval([datatype '(nidat.img)']);
                end

                sandbox.hdr = nidat.hdr;
                scaninfo.hdr = nidat.hdr;
                sandbox.hemodynamics = [];
                sandbox.onset = [];
                sandbox.intrial = [];
                clear nidat;
                waitstr = (['reading file 1 of ' num2str(filnum)]);

                for x = 2:filnum
                    waitbar(x/filnum,h,waitstr)
                    %nidat = avw_img_read_mt(strcat(file_names{x}));
                    nidat = load_nii(strcat(file_names{x}));
                    waitstr = (['reading file ' num2str(x) ' of ' num2str(filnum)]);
                    disp (['Loading file ' file_names{x}]);


                    %nidat = eval([datatype '(nidat.img)']);
                    %idat(:,:,:,x) = nidat;
                    idat(:,:,:,x) = nidat.img;
                    clear nidat;
                end

            case '2dseq'
                % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Load Bruker 2dseq images
                % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                waitstr = ['Loading ' num2str(filnum) ' files'];
                h = waitbar(0,waitstr);
                [idat sandbox scaninfo] = fsb_load_2dseq(pathstr,file_names);

                sandbox.name = '2dseq file';
                waitstr = (['reading file 1 of ' num2str(filnum)]);

            case 'SPM'
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Load files indexed in a SPM struct
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                disp ('Loading SPM.mat');
                tempdata =load(file_names{1});
                SPM = tempdata.SPM;
                disp ('Loaded SPM.mat');
                file_names_spm = SPM.xY.P(:,1:end-2);
                filnum2 = size (file_names_spm,1);
                waitstr = ['Loading ' num2str(filnum2) ' files'];
                h = waitbar(0,waitstr);
                nidat = avw_img_read_mt(strcat(file_names_spm(1,:)));            
                datatype = sub_get_data_range(nidat.img, datamode);
                idat = eval([datatype '(nidat.img)']);
                clear nidat;

                for x = 2:filnum2
                    waitbar(x/filnum2)
                    nidat = avw_img_read_mt(strcat(file_names_spm(x,:)));

                    nidat = eval([datatype '(nidat.img)']);
                    idat = cat(4,idat, nidat) ;
                    clear nidat;

                end

                if isfield(tempdata.SPM,'sandbox');
                    sandbox = tempdata.SPM.sandbox;
                else
                    sandbox = fsb_import_spm(SPM);
                    sandbox.disp_onset =1;
                end

                disp([ num2str(filnum) ' loaded in ' num2str(t) ' seconds' ]);

            case 'mat'
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Load sandbox/idat file pair
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                tempdat = load(file_names{1});
                disp (['loading file ' file_names{1}]);

                for ii = 1:filnum;
                    s = whos ('-file', file_names{ii});
                    temp(:,ii) = s.size;
                    [pathstr{ii}, name{ii}, ext, versn] = fileparts(file_names{ii});
                end
                
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Try to load Monline tcImg format
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                
                if isfield(tempdat,'tcImg')
                    [idat sandbox] = fsb_import_monline(file_names);
                end
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Check for idat field
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                if isfield(tempdat,'idat')

                    temp2 = temp(1:3,1);

                    if length(temp)>3;
                        temp2(4) = sum(temp(4,:));
                    else
                        temp2(4) = 0;
                    end

                    clear idat idat_spare imean;

                    [a b c d] = size(tempdat.idat);
                    if d>5000;
                        for x =1:size(file_names,2)
                            tempdat = load(file_names{x});
                            [aa(x) bb(x) cc(x) dd(x)] = size(tempdat.idat);
                        end
                        if [max(aa) max(bb) max(cc)]-[min(aa) min(bb) min(cc)]>0;
                            misfit = 1;
                        else
                            misfit = 0;
                        end
                    else
                        misfit = 0;
                    end

                    disp ('Loading image information...');
                    tempdat = load(file_names{1});
                    datatype = sub_get_data_range(tempdat.idat,datamode);
                    tempdat.idat = eval([datatype '(tempdat.idat)']);

                    curlength = size(tempdat.idat,4);
                    if size(file_names,2)>1 && curlength>1;
                        
                        idat = zeros(temp2(1),temp2(2),temp2(3), temp2(4),datatype);

                        if misfit ==1;
                            idat = fsb_pad_idat(tempdat.idat,0,[max(aa) max(bb) max(cc)]);
                        else
                            idat(:,:,:,1:temp(4,1)) = tempdat.idat;
                        end
                    else
                        idat = tempdat.idat;
                    end

                elseif isfield(tempdat,'tcImg')
                    idat = fsb_read_tcimage(file_names);
                    sandbox.name = 'multiple files';
                    scaninfo = [];
                    sandbox.hdr.dime.pixdim(2) = 99;
                    sandbox.hdr.dime.pixdim(3) = 99;
                    sandbox.hdr.dime.pixdim(4) = 99;
                    return
                else
                    idat = [];
                end

                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Check for sandbox field
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                if isfield(tempdat,'sandbox');

                    tempdat = tempdat.sandbox;
                    sandbox = tempdat;

                    try
                        tempdat = fsb_volmat(tempdat);
                    catch
                        disp ('Nothing new here')
                    end
                    try
                        sandbox.sb{1}= tempdat;
                        sandbox.subname{1}= name{1};
                        sandbox.time{1} = sandbox.timemat;
                        sandbox.time{1}.volmat = tempdat.volmat;
                        sandbox.volmat = tempdat.volmat;
                        sandbox.time{1}.curlength = length(sandbox.hemodynamics);
                    catch
                        disp('sandbox struct not complete');
                    end

                    if size(sandbox.intrial,2)< 4 ;
                        sandbox.intrial(:,4)= 1; % Run number
                    end

                    if size(sandbox.intrial,2) < 5;
                        sandbox.intrial(:,5) = 1; % Session number
                    end

                    if size(sandbox.intrial,2) < 6 ;
                        sandbox.intrial(:,6)=1; % Subject number
                    end

                    if size(sandbox.intrial,2)< 7
                        sandbox.intrial(:,7) = 1:size(sandbox.intrial(:,2));
                    end
                    disp ('Sandbox struct found, loading sandbox variables...');
                    [sandbox prms] = fsb_import_sandbox(sandbox,idat);

                else

                    sandbox.hdr.dime.pixdim(2) = 99;
                    sandbox.hdr.dime.pixdim(3) = 99;
                    sandbox.hdr.dime.pixdim(4) = 99;

                end

                if ~isfield(tempdat,'scaninfo');
                    scaninfo = [];
                end


                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                % Load more files
                %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                h = waitbar(0,waitstr);

                for x = 2 : size(file_names,2);
                    waitbar(x/size(file_names,2));
                    tempdat = load(file_names{x});
                    disp (['loading file ' file_names{x}]);

                    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                    % Try to load more idat files
                    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                    if isfield(tempdat,'idat')

                        if misfit ==1;
                            tempdat.idat = eval([datatype '(tempdat.idat)']);
                            p_idat = fsb_pad_idat(tempdat.idat,0,[max(aa) max(bb) max(cc)]);
                            idat = cat(4,idat, p_idat);
                        else
                            if size(temp,1) == 4; % if whole functional datasets are concatenated
                                idat(:,:,:,curlength+1:curlength+temp(4,x))= tempdat.idat;
                                curlength = curlength+temp(4,x);
                            else % if only 3D images are concatenated
                                idat(:,:,:,curlength+1)= tempdat.idat;
                                curlength = curlength+1;
                            end


                            tempdat.idat =[];
                        end
                    end

                    %~~~~~~~~~~~~~~~~~~~~
                    % Try to load more sandbox files
                    %~~~~~~~~~~~~~~~~~~~~
                    if isfield(tempdat,'sandbox');
                        try
                            tempdat.sandbox = fsb_volmat(tempdat.sandbox);
                        catch
                            disp ('Sandbox struct already in new format')
                        end

                        try
                            tempdat = tempdat.sandbox;
                            sandbox.sb{x}= tempdat;
                            sandbox.subname{x}= name{x};
                            sandbox.time{x} = tempdat.timemat;
                            sandbox.time{x}.curlength = curlength;
                            sandbox.time{x}.volmat = tempdat.volmat;

                            onsets = sandbox.volmat.evt_onsets;
                            onsets2 = tempdat.volmat.evt_onsets;
                            durs = sandbox.volmat.evt_dur;
                            durs2 = tempdat.volmat.evt_dur;

                            for ii = 1:size(onsets,2)
                                onsets_o = onsets{1,ii};
                                onsets_n = onsets2{1,ii};
                                dur_o = durs{1,ii};
                                dur_n = durs2{1,ii};
                                dur_n = [dur_o;dur_n];
                                onsets_n = onsets_n + sandbox.time{x-1}.curlength;
                                onsets_n = [onsets_o;onsets_n];
                                onsets_n = sort(onsets_n,'ascend');
                                onsets{1,ii} = onsets_n;
                                durs{1,ii} = dur_n;
                            end

                            sandbox.volmat.evt_dur = durs;
                            sandbox.volmat.evt_onsets = onsets;
                            disp ('Sandbox struct found, loading sandbox variables...');
                        catch
                            disp('Sandbox struct cannot be loaded');
                        end
                    end

                    %~~~~~~~~~~~~~~~~~~~~
                    % Try to load hemodynamics files
                    %~~~~~~~~~~~~~~~~~~~~
                    try
                        sandbox.hemodynamics = [sandbox.hemodynamics;tempdat.hemodynamics(:,1:size(sandbox.hemodynamics,2))];
                        sandbox.disp_onset = 1;
                        sandbox.hrf_pred=1;
                        disp ('Loading hemodynamics data...');
                    catch
                        disp ('no hemodynamics data loaded');
                        sandbox.hemodynamics = zeros(1);
                        sandbox.disp_onset = 1;
                    end

                    %~~~~~~~~~~~~~~~~~~~~
                    % Try to load hemodynamics_deriv files
                    %~~~~~~~~~~~~~~~~~~~~
                    try
                        sandbox.hemodynamics_deriv = [hemodynamics_deriv;tempdat.hemodynamics_deriv];
                        disp( 'Loading hemodynamics_deriv data...');
                    catch
                        disp( 'Error loading hemodynamics derivatives');
                    end

                    %~~~~~~~~~~~~~~~~~~~~
                    % Try to load onset files
                    %~~~~~~~~~~~~~~~~~~~~
                    try
                        sandbox.onset = [sandbox.onset;tempdat.onset];
                        sandbox.disp_onset = 1;
                        sandbox.hrf_pred=1;
                        disp( 'Loading onset file...');
                    catch
                        disp( 'Error loading onset file');
                    end
                    
                   %~~~~~~~~~~~~~~~~~~~~
                    % Try to load sensdata files
                    %~~~~~~~~~~~~~~~~~~~~
                    try
                        sandbox.sensdata = [sandbox.sensdata;tempdat.sensdata];
                        disp( 'Loading sensor data information...');
                    catch
                        disp( 'Error loading sensor data information');
                    end

                    %~~~~~~~~~~~~~~~~~~~~
                    % Load intrial files
                    %~~~~~~~~~~~~~~~~~~~~
                    try
                        if size(tempdat.intrial,2)<3 && ~isempty(tempdat.intrial);
                            sandbox.intrial = [sandbox.intrial; [tempdat.intrial(:,1) (ones(size(tempdat.intrial,1),2))*x]];
                            disp( 'Loading intrial file...');
                            intrial2 = sandbox.intrial(:,1);
                            intrial2(intrial2==0.01)=0;
                            sandbox.intrial(:,3) = bwlabeln(intrial2(:,1));
                        else
                            tempdat.intrial(tempdat.intrial(:,3)==0,3)=(-max(sandbox.intrial(:,3)));
                            tempdat.intrial(:,3) = tempdat.intrial(:,3)+max(sandbox.intrial(:,3));

                            if strcmp(name{x}(1:6),name{x-1}(1:6))==0
                                keepname = keepname+1;
                            end

                            if strcmp(name{x}(1:3),name{x-1}(1:3))==0
                                subjectname = subjectname+1;
                            end

                            if size(tempdat.intrial,2)<4;
                                tempdat.intrial(:,4)= x; % Run number
                            end

                            tempdat.intrial(:,5)=keepname; % Session number

                            if size(tempdat.intrial,2)<6
                                tempdat.intrial(:,6)=subjectname; % Subject number
                            end

                            if size(tempdat.intrial,2)<7
                                tempdat.intrial(:,7) = 1:size(tempdat.intrial(:,2));
                            end

                            sandbox.intrial = [sandbox.intrial; tempdat.intrial];

                        end

                        disp (['Number of correct trials: ' num2str(max(sandbox.intrial(:,3)))]);
                    catch
                        disp('Error loading intrial files')
                    end

                    clear tempdat;
                    waitbar(x/filnum);
                end
        end
    end

    close(h);
    sandbox.name = [name{1} ext];


end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Use this function to determine the datatype of the image data to load
% Data with a variance lower than 1000 or with values lower than zero
% will automatically be converted to double
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function datatype = sub_get_data_range(image_data,datamode)
        data_min = min(image_data(:));
        data_max = max(image_data(:));
        datatype = 'int16';
        if data_max - data_min < 1000
            datatype = 'double';
        end
        if data_min < 0
            datatype = 'double';
        end
        if data_max > 32768
            datatype = 'double'
        end
        
        if datamode == 0
            datatype = 'int16';
        end
    end

end

